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ABSTRACT 


A computer program that calculates absorption spectra for slant paths 
through the atmosphere is described. The program uses an efficient con- 
volution technique (Romberg integration) to simulate instrument resolution 
effects. A brief information analysis is performed on a set of calculated 
spectra to illustrate how such techniques may be used to explore the qual- 
ity of the information in a spectrum. 





I. INTRODUCTION 


Radiation which has passed through and been modified by the atmo- 
sphere contains a wealth of information about atmospheric conditions. 

This includes information about the temperature and density structure of 
the atmosphere as well as the abundance and distribution of constituent 
gases and aerosols. In order to retrieve this information, however, an 
accurate model is needed which describes the dependence of absorption 
spectra on the atmospheric conditions and on the instrument characteris- 
tics. The first dependence has received considerable attention, both 
theoretical and computational, but the second has received much less. 
Equivalent width techniques and transmittance averaging are common ways 
of extracting information from spectra when instrument effects are un- 
known or poorly understood. However, because these effects mask under- 
lying atmospheric effects, knowledge of the instrument behavior, when 
properly included in the data analysis, allows extraction of more detailed 
information than can be obtained in the absence of such knowledge. 

Nonlinear least- squares 1 is one data analysis technique which can 
include all relevant influences on the data. It involves a great deal of 
computation, however, so efficient computer codes must be developed in 
order to take advantage of the least- squares techniques. This report de- 
scribes a program for calculating absorption spectra that was specifically 
written to efficiently include instrument effects. It is therefore suit- 
able for use in least-squares analysis. 

In constructing a model for any experiment it is important to know 
what type of information the experimental data are likely to contain. 

Such knowledge allows the level of sophistication for the model to be 
balanced between computational complexity and interpretative accuracy. 

The third chapter of this report addresses this question and illustrates 
a variety of simple techniques for information analysis with a set of at- 
mospheric absorption spectra. 
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II. FUNDAMENTALS OF ATMOSPHERIC TRANSMITTANCE 


A beam of radiation propagating through the atmosphere is modified 
through the action of many different processes. Among these are scatter- 
ing, reflection, refraction, polarization changes, atmospheric emission, 
and absorption. This report is only concerned with the modifications due 
to refraction and absorption by gas molecules in the atmosphere, and, in 
particular, that portion of the absorption due to rotational and vibra- 
tional transitions of the molecules. Each transition has associated with 
it a frequency position, width, intensity and profile. The Voigt pro- 
file 2 ’ 3 ’^ 3 is assumed throughout this report, although other shapes can 
be incorporated if desired. 

A. HOMOGENEOUS PATHS 

If radiation with a spectral irradiance E(v) is incident on a homo- 
geneous gas sample, the spectral radiant exitance M(v) is given by 4 

M(v) = E(v) • t(v) , (1) 

where v is the frequency of the radiation measured in cm -1 . The modifying 
factor, t ( v ) , in Eq. (l) is called the spectral transmittance, which is 
referred to simply as transmittance in this report, it being understood 

that the transmittance is a function of frequency. 

: 

For a narrow homogeneous sample of thickness di, Bouguer's law 5 
states that the radiation lost in the sample is linear in E(v) and d & . 
Beer's law 5 states that the radiation lost is also linear in the concen- 
tration N of the absorber, so that 

dE(v,i) = -k(v) E(v,i) Ndi , (2) 

where k(v) is called the absorption coefficient. When N is measured in 
molecules/cm 3 and di! in cm, k(v) is measured in cm 2 /molecule . When N is 
independent of path length, Eq. (2) can be integrated over the entire path 
length, L, through the sample to give 

M(v) = E(v,L) = E(v,0) exp[-k(v) NL] . (3) 

Consequently, the transmittance is given by 

t(v) = exp[-k(v) NL] . (4) 

Equation (2) is based on independent absorption by each absorbing molecule. 
It does not apply to self -broadening absorption, therefore, in which k can 
also depend on N. The factor in the exponential of Eq. (4) is called the 
absorbance. 

If there are M different absorbers with concentrations % , Ns, ..., Njq 
independently absorbing in the sample, the total absorbance A(v) is 
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( 5 ) 


: M . 

A(v) = E iysm(v) L 

m=l 

where l%(v) is the absorption coefficient for the mth absorber. As stated 
above, this report assumes that the absorption is due to many individual 
absorption lines, each with absorption coefficient kmj(v). Consequently 
Eq. (5) becomes , 

; i 

' M ' Jm 

A(v) = L • XI N m ‘ E > (5 d ) 

i m=l j J=1 

I s 

where J m is the total number of lines for the mth absorber. 

i ! ! 

B. ATMOSPHERIC PATHS 

! 

| Radiation paths in the atmosphere are usually not homogeneous. For 
such paths, Eq. (2) can still be integrated to give E as a function of , 
path length, provided the variation of k and N with path length is known'. 
The analog of Eq. (5) for atmospheric paths is I 

j : 

I rL M j Jm 

i A(v) = / X] %i(i) E ^mj( v >-0) ^ • (6) 

j J 0 m=l j=l 

i 

It is convenient to define the mixing ratio a m (T) for the mth absorber as 
^ am(jg) = %(£)/N(i) , 

t , ; 

where N(i) is the mean density of the atmosphere as a function of path | 

length. | 

| Defining a new variable u(i), the so-called airmass along the path, by 
j du = N(£) d£ , I 

f ' I 

;< . i 

Eq. (6) becomes j | 

i; ’ 

jj HJ M ! J m ' 

A( v) = / E %i( u ) E k mj( v J u ) du » (7) j 

J J 0 m=l 1 j=l j 

where I 

H i ! 

u = jf^NU) dl 

j j 

is the total airmass of the path. J ; 
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C. REFRACTION EFFECTS 


This report uses a program written by Treve 7 ’ 8 and Snider 9 to deter- 
mine the path of radiation through the atmosphere. The necessary geometry 
is shown in Fig. 1, assuming that the sun is the radiation source. The 
observer is at a height above the surface of the earth at latitude <t> . 
The direction of propagation of the radiation makes ah angle 6 1 with the 
observer's local zenith. A straight line drawn from the observer to the 
source makes an angle 9 with the local zenith and has a minimum height 
HjJiin above the surface of the earth. Because of refraction, the true 
radiation path is not straight and has a different minimum height %i n . 

The actual path taken by the radiation depends on the density profile of 
the atmosphere as a function of height. This profile is calculated by 
the program from an imput temperature profile and base level pressure by 
integrating the hydrostatic equation and assuming an ideal gas law for the 
relationship between pressure, temperature, and density. 10 The input 
profile is in the form of many thin. atmospheric layers. The program uses 
Snell's law and the Lorentz-Lorenz relation to perform a ray- tracing 
through the layers. The output consists of a series of airmass values, 
effective temperatures, and pressures, one for each layer crossed by the 
ray. 


In order to use these to calculate the spectrum, the integral in 
Eq. (7) is approximated by a sum according to 


/ : U M Jm 

Y am(u) Y k mj( v 5 u ) 

0 m=l j=l 

IL M • J m 

- Y Y Y L a mkmo( ,v )J s (8) 

i=l m=l j=l 1 

where 1^ is the number of terms in the sum -which replaces the integral in 
El. (7). 

The form of [a m ^mj( v )]. 1S determined by assuming that the product 

a m (u) k ffl j(v,u) for fixed' v can be written as a function of the pressure 
P(u) and temperature T(u) along the path. We therefore write 


/ %(u) l%nj(v,u) du = / f[P(u),T(u)] du . 
Expanding f in a Taylor series about some value, P 0 , To, 


f(P,T) = f(Po,To) + || 


(P-Po) 


Po,T 0 


(T - To) + {second-order terms) 
Po,To 


Substituting this expression into the integral 
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6 = "true" or astronomical zenith angle 
0 ' = apparent zenith angle 
H 0 b = observer height 
H^in = unrefracted minimum height 
Hjnin = true minimum height 
R(0) = radius of earth at latitude 0 


Fig. 1 - Observer- source geometry 




/ f[P(u) jT(u) ] du = f(Po,To) / du 

+ ||| f (P - Po) du + / [T(u) - To] <** 

Po j To Po, To 


where 


+ (second-order terms]' 
[f(Po,T 0 ) + 


df /p p \ , ’^f 

.5?-- (P - Po ? + ,§T 


(T - To) 


\ • / 


Po , To i. Po , To 

"+■ .{second-order terms) ]/ 


du* 


P = / P(u) du/ / du 


and 


(9) 


T = / T(u) du// du . 

But the expression in the brackets, to first-order, is just f(P,T) . 
Consequently, provided the variation of P(u) and T(u) is sufficiently small 
that the first-order terms in the Taylor series are adequate, we may write 

* 

/ f[P(u),T(u)] du - f (P,T) Au , (10) 

where T and P are defined as in Eq. (9) . 

If the atmospheric layers are sufficiently thin, we are therefore 
justified in replacing the integral- in Eq. (7) with the sum in Eq. (8) 
provided 

[a m k mJ (v)]. = a m (Pi,Ti) k raj (P t ,f t , v) , (11) 

where 

Pi = f P(u) du/AUj 
-layer • i 

and 1 (12) 

T*i = J T(u) du/AUi . 
layer i 
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Equations (11) and (12) are close to the familiar scaling approxima- 
tions . 11 The same argument applies to any other variables that effect the 
a • k product ( such as mixing ratio dependence in k) . One can therefore 
approximate the inhomogeneous atmosphere with a series of homogeneous 
layers whose properties are given by averages over the true inhomogeneous 
variation. 


D. LINE SHAPE 


As stated above, this report assumes that the absorption coefficient 
k m< 3 (v) follows the Voigt profile. This profile combines the effects of 
collision broadening (which is dominant in the troposphere) and Doppler 
broadening (which is dominant in the upper stratosphere) . It is given by 


km j ( v 5 Lj l) — ^ 


r°° exp(-t a ) dt 
^ + (x - t ) 2 5 


(13) 


where 

komj = V^’2/« > - ' 

a D 

y = — - V^2' , and 
a D 

X = — V ° J ” yj^n2 . 

a D 


The line intensity Sj m is defined as Sj m = /q kj m (v) dv . 
From Quantum - Mechanics 12 this - quantity is given by> 


Sj m 


8 rt 3 v 0i m g"(l - e~ hcv Oj'”/ kT ) 
3 hc QjT) 


exp[(-E„j )/kT] • 


(14) 


In this expression, v 0 j m is the line center frequency for the mjth 
line with lower state energy E^j and upperstate energy E^j = hcv 0m j + E"j 
g" is the degeneracy of the lower state and Q m (T) is the partition funcl 

tion. In the infrared, -the expression eX ^ can usually be 

approximated by l/Q ro t( T ) where Q ro t(T) , the rotational partition function 
is proportional 1 3 to T for linear molecules such as CO, N 2 O, CO 2 , and NO. 
It is proportional to t 3//2 for molecules like O 3 , H 2 O, and CH 4 . The term 
is the square of the dipole matrix element. 
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The temperature dependence of S jra is made explicit by writing it as 


m 


= Soj 


J ra 


(^f exp[-E«j(T 0 - T)/0.6946TT o ] 


(15) 


where p = 1 or 3/2 depending on the molecular structure. The above expres- 
sion assumes that E is measured in cm -1 and T in Kelvin. S 0 j m is the value 
of Sj m at some standard temperature, To. 


The quantities and cc-p in Eq.' (13) are the Lorentz width and Dop- 
pler width of the line, respectively. These are given by 


a D 


vom.i l2kT 
e \ M n -*- 


and 

«L = Qbmo (cf) ■ 


(16) 


M m is the mass of molecular species m; ao m j is the Lorentz width at 
some standard temperature and pressure. To and Po; and 7 is an empirical 
term describing the temperature dependence of o^. This report makes the 
standard assumption that 7 = l/2 although other values are often more 
accurate for bands of specific molecules. The necessary line parameters 
were taken from the AFGL line listing. 14 

E. INSTRUMENT EFFECTS 

When Eq. (13) and (8) are combined, they constitute a mathematical 
model for the frequency make-up of a beam of radiation as it propagates 1 
through the atmosphere. Since real spectrometers have finite resolution, 
the signal observed at frequency setting v is a weighted sum of the radia- 
tion incident on the instrument. This is described by writing 

T con( v ) = f'o CT ( V > V ’) T mon( v ’) dv ’ > 

where T con (v) is the observed transmittance, % 0 n( v ') 1S 6he m° noc i^ oma li c 
transmittance, and cr(v,v') is the instrument spectral response function 
(or "slit" function) with 

Iq v( v >V) dv' =1 . 

It is often assumed that cr can be written as a function of (v - v ' ) , 
and many different forms are used. For grating instruments, triangular 
and Gaussian shapes are common, 15 and for interferometers sine, sine 2 , and 
various Bessel functions have been tried, 16 depending upon the apodization. 

When calculated and observed spectra are to be compared visually, the 
exact form of a is often unimportant. However, when more sensitive 
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techniques, such as least-squares analysis, are used, instrument effects 
must be modeled with a precision at least better than the precision of 
the observed spectral values. This imposes an additional set of computa- 
tional restrictions on spectrum calculation routines, especially when the 
convolution integral in Eq. (17) must be evaluated numerically, as is 
usually the case. While considerable effort has gone into devising ef- 
ficient algorithms for calculating monochromatic transmittances, little 
work seems to have been done on the special problems associated with cal- 
culating convolved transmittances. 

Since one of the problems this report seeks to address is the appli- 
cability of least- squares data analysis to atmospheric transmission spec- 
tra, the precision of convolved transmittances has been a central concern 
from the start. To see how this concern dictates the overall calculation 
strategy, consider how Eq. (17) is evaluated in practice. Because of the 
complex frequency dependence of -r(v'), the integral must be evaluated 
numerically; i.e., replaced by a weighted sum. Each term in the sum re- 
quires a value of T mon and one of cr. If N c is the number of terms in the 
sum and N s is the number of observed spectral values, then N c • N s mono- 
chromatic transmittance values are needed to calculate the convolved 
transmittance, -r con . 

A typical frequency position may have contributions from 20 different 
lines and require 20 atmospheric .layers to model the atmosphere (Mankin 17 
used ~1000 lines over a 10 cm -1 region; Snider and Goldman 10 used 197 
layers) . In addition, if the resolution is &vp> a typical frequency spacing 
for observed spectra is Svp/2, Consequently for a 200 cm _1 -long spectrum 
at moderate resolution, 6vp =0.10 cm -1 (some interferometers have 10 
times better resolution) , a total, of 4000 observed points must be calcu- 
lated. Assuming that the most time-consuming part of the calculation is 
the evaluation of the Voigt profile, a total of N c X 20 X 20 X 4000 = 

1.6 x 10 6 U C Voigt evaluations must be performed. Running time is dependent 
on machine characteristics, but Pierluissi and Vanderwood 18 quote a time 
of ~3 seconds for 100 evaluations on an IBM 360/65. This gives 
4.8 X 10 4 • N c seconds of CPU time in the Voigt routine alone. Clearly 
N c must be made as small as possible. 

1. Convolution Strategies 


The numerical evaluation of the integral in Eq. (17) is represented 
schematically in Fig. 2. To each observed frequency position Vj there 
corresponds a set of frequency positions {v'}i at which monochromatic 
transmittances must be calculated. Clearly the number of points v' can 
be reduced by making some of the { v ' ) i +1 points overlap some of the ,{v')i 
points. However, because the number of { v ' } ± values needed for a given 
level of accuracy e in T con (vi) depends upon the form of the monochromatic 
transmittance, the density of {v'Ji points will be a function of for a 
given e . The exact density depends on the numerical integration rule em- 
ployed, but usually this takes the form of bounds on the magnitude of high- 
order derivatives of T mon over the convolution frequency range. The use 
of an iterative quadrature routine can often eliminate the need for 
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monochromatic 



convolved 


Fig. 2 - Schematic representation of Eq_. (17) 



evaluating these derivatives. Also, if the T con ( Vj ) values are calculated 
in order of increasing [or decreasing] frequency position, values of T mon 
at frequencies lower [higher] than the ith range of v' values need not be 
saved. This reduces storage requironents . 

Of the many available iterative routines, am equally spaced abscissa 
Romberg integration algorithm was selected. Gaussian rules require fewer 
[v 1 ] points for a given e, but are not well adapted to iterative evalua- 
tion. Neither are they well suited for overlapping (v'} positions for 
successive v £ values. Adaptive Simpson's rule techniques are good at 
varying the v' -density to meet the accuracy requirements and could be writ- 
ten so as to use overlapping £ v ' ) positions. Unfortunately they are much 
more complicated than Romberg schemes, which proved to be complicated 
enough. Better programmers" than the_ author of this report- may wish to use 
them, however. 

2. Romberg Integration 


To use classical Romberg integration to evaluate the integral 

I = f b f (X)dx 

J a 


the interval [a, b] is divided into 2^° equal pieces. The function f is 
evaluated at each of the 2^° + 1 points, 

Xi = a + (i - l)(b - a)/2 k o i = 1, . . . , + 1 , 

and the sum, 

T e. jSo) = [ 1/2 fUi : a) + f(x 2 ) ♦ . . . 

+ f(x 2 ko) + V 2 f ( x 2 k o,+l = b )] ’ ( l8 ) 


forms the first element of a table of approximations to I called the 
Romberg T- table. Equation (18) .is easily recognized as the trapezoidal 

" C' r - ' (kn) - . - - ~ - 

rule. Next, .ko is replaced by k L = ko- + 1 and To is calculated from 
Eq. (18) with k 0 replaced by k L . This value forms a second estimate for 

I and is written next to T^^°^ in the top-most row of the table, as shown 

- - " ~ ■ / V"\ _ . ‘ ' 

in Fig.- 3. From these two estimates'/a third estimate, T^ , is formed 
according to 



( 19 ) 
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rp(ko) T (kl) mOfe) . . 

o o o 


rn(ko) rp'(kl) 

1 1 


rjlCko) . . 

2 


Fig. 3 - The Romberg T-table 
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are then compared, and if 


The values of T 


(ko) 


and T 


(ko) 

1 


rp(ko) m(ko) 

i l o 

the process stops and outputs as the final approximation to I. If 

the convergence criterion in Eq. (20) is not satisfied, the table construc- 
tion continues with ^ , Tq^ 3 '* and so on until two successive elements 
along the diagonal of the table do satisfy the convergence criterion. 


< e 


( 20 ) 


Davis and Rabinowitz 19 discuss a modification of the classical Romberg 
scheme which requires calculation of a second table, the Romberg M- table. 

In this scheme, after Tq^°^ is calculated, the quantity 


M. 


(ko) 


2ko 

) - a v* 

2ko £l 


'[a + - |)(b - a)/2 ko ] 


( 21 ) 


is calculated. 
The values of T 


Equation (21) is easily recognized as the midpoint rule. 


(ko) 


and M, 


(ko) 


are then compared and if 


(ko) _ (ko) 
o o 


< 


( 22 ) 


the process stops and outputs 

= I ( T o ko) + (23) 

as the final approximation to I. If the new convergence criterion [Eq. 
(22)] is not met, the next diagonal element of the T-table is calculated 
from 


T (k) 

m 


M 


(k) 

‘m-l 


(* • ^ - tl) 


k m - 1 


(24) 


The next element of the M-table, M^ 1 ^ , is calculated from Eq. (21) ' 
with ko replaced by ki = ko + 1. The two M values are then combined as 

(kn) 

in Eq. (19) to form ° . The convergence criterion is then applied to 
and T^ 1 ^ . If the criterion is satisfied, these two are combined 


as in Eq. (23) and output. Otherwise the process continues halving the 
interval [a,b] and recursively generating the T- and M- tables until two 
corresponding diagonal elements satisfy the convergence criterion. 


Those referring to Davis and Rabinowitz 19 should note that Eq. (6.3.9) 
in their book is incorrect in some editions, as can be verified from their 
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Eq. (6.3.8) and the equation below (6.3.2). Also, the FORTRAN version of 
this technique, incorporated into the program in Appendix A, has a maxim um 
number of interval halvings in the search for convergence. When this 
l imi t is exceeded, the output is taken from Eq. (24) rather than Eq. (23). 

Incorporating the Romberg Table evaluation in the spectrum calcula- 
tion routine required less than 15 FORTRAN statements. The logic needed 
to save T mon values required nearly 10 times as many statements. The de- 
tails are given in Appendix A. 

F. VERIFYING CALCULATION BEHAVIOR 

The program which evaluates Eq. (17) , described in Appendix A, is 
composed of three basic parts: (l) evaluation of the atmospheric path 

(subroutine SNIDER), (2) setting up approximate layer parameters (subrou- 
tine PATHST) , and (3) spectrum calculation (subroutine NSPEC) . These were 
first checked individually and then as a whole. 

'A^speci f i c_.ob server^ s our c e geometry was "selected. __This._had’ an. observer 
at 45°;- N latitude, jij+!^“km height,, .and an. .observed' ,zbnithjangle (0' in' Fig. 1‘) 
of 95.2°. The source was assumed to be above the atmosphere. Unfortu r 
nately, the atmospheric model used by Snider and Goldman 10 could not be 
obtained so the January reference atmosphere of Cole and Kantor 20 was used 
instead. The program's analysis of the path is shown in Table 1’ along 
with the associated results taken from Ref. (10). Also shown are the re- 
sults using a 190-layer atmosphere taken from the U.S. Standard Atmosphere , 
1976 . 21 The agreement shown in the table and the proper behavior of the 
subroutine to changes in zenith angle, observer height, and ground-level 
pressure indicated that the subroutine was functioning properly. 

The output for the 190-layer model included 198 sets of (ui, Ti, Pi) 
values; 59 values for layers crossed in the portion of the path going from 
the observer in layer #99 to the minimum height in layer #4l, and 139 
values for layers crossed in the portion from the minimum height to the 
top of the atmosphere. The average layer thickness in this model is 
~.4 km. Examination of these 198 sets showed that of the 1.200099 x 10 26 
mol/cm 2 in the path, 3 . 13^1 X 10 25 , or 26$, were in the 0.35 km thick por- 
tion of layer 4l that contains the lowest segment of the path. This dem- 
onstrates the height selective nature of the limb geometry.. 

Compare this to the result obtained by Schmidt 22 that a uniform error 
of. 0.5 km .in the minimum height of a. limb emission spectrum produced a : 

21$ rms error in retrieved mixing ratios (using simulated balloon data) . 

The next highest layer, #42, contained 1.4154 X 10 25 mol/cm 2 , or 
11.8$. A full 60$ of the air mass lay between 18.45 and 20.49 km. 

Next the Voigt function evaluation routine, due to Hui, Armstrong, 
and Wray, 23 was checked against the tables of Young. 24 It gave good 
agreement over the tabulated range of x and y values. A set of sample 
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Table 1. Comparison of path characteristics. 



Ref. (20) 

Ref. (10) 

Ref. (21) 

Apparent Zenith angle 
(degrees) 

- • - 95.. 2 

95ii2' ' 

95 *.2 ' 

Lower height (m) 

45000.0 

45000.0 

45000.0 

Upper height (m) 

89999.9 

99000.0 

82000.0 

Minimum height (m) 

18454.5 

18434.1 

18450.0 

Chapman airmass 

6.168 

7.098 

6.385 

Optical airmass 

5.422 

6.085 

5.573 

HDIFF (m) 

1000.5 

1138.5 

1030.0 

Refraction (degrees) 

0.11255 

0.12788 

0.11579 

Astronomical Zenith 
angle (degrees) 

95.31 

'95.33 

95.32 

P(EFF) (w/m 2 ) 

4518.41 

5494.93 

4999.92 

t(eff) (k) 

216.26 

219.54 

218.93 

•s ^ 

P(TM) (N/rn 2 ) 

7890.17 

9066.33 

‘ 

. 
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line parameters were then chosen. These are shown in Table 2. A uniform 
mixing ratio of 0.1 ppm was used to calculate the absorbance in each 
layer . 

The total absorbance was 1.231378 at v = v 0 = 2090. 796 cm -1 , cor- 
responding to a transmittance (no convolution) of 0.29189. The absorbance 
in layer 4l was 0.197407, or ~16$, and the absorbance between 18.45 and 
20.49 km was 0.456708, or 37$. These values show that the absorbance at 
line center is less peaked with height than the airmass. For Lorentz-like 
lines, the absorbance at v Q is Su/cmt. Consequently the decrease in u with 
height is offset by the decrease in a (because of decreasing pressure) and 
the absorbance is "smeared out" with height. The degree of Lorentz-like 
behavior for Voigt lines is determined by y in Eq. ( 13 ) . Large y implies 
a high degree of Lorentz-like behavior; small y implies mainly Doppler- 
like behavior. The Lor entz- Doppler crossover point, y =* 1, occurs at layer 
63 for the sample line and atmospheric profile used; consequently, the 
observed "smearing" is reasonable. 

The above observations suggest that high-resolution spectra should 
have better height sensitivity than low-resolution spectra because the 
same effect which "smears" the line center absorbance sharpens the wing 
absorbance, which goes as Sua. In low-resolution spectra the observed 
signal is a mix of the two and has intermediate height sensitivity. High- 
resolution spectra are better to use to study the height- sensitive wing 
absorbance. 

A second consequence of the "smearing" is found when the fine layers 
are combined into thicker layers. Let consecutive layers __ 
ii, £2, . . . , Ji n with parameters (u^Ti^), . . . , (u n ,T„,P„) be com- 
bined into one layer with parameters 

n 

n = £ J 
i=l 

T = £ U i T i/ U > 311(1 (25) 

P = £ UiPi/u • 

The new layer represents a less accurate approximation to the atmosphere 
than the original layers. In order to save computation time, it is 
desirable to use the minimum number of layers to calculate the spectrum. 
Subroutine SNIDER is sufficiently fast, however, that no such concerns 
apply to it. To test the sensitivity of the calculated spectrum to the 
layer coarseness, the original 139 sets of (u,T,P), one for each layer 
crossed, were combined into segments of approximately equal airmass . Each 
new layer had an airmass as close to, but less than, U/Nge as possible, 
where Nge is some integer. The calculated absorbances at three 'different 
frequencies for several different values of Nge are shown in Table 3. It 
is clear from these results that more layers are required at line center 
than in the wings. This is because the result of decreasing Nge is to 
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Table 2. Sample line parameters* 


v 0 = 2090.796 cm 1 
S 0 = 0.189 X10“ 2 ° mol -1 cm 
Oo = 0.110 cm -1 
E" = 272.220 cm -1 
m = 48 amu 


Corresponding to the (6,13)-(6 } l4) line of the v x + v 3 band of 16 0 3 . 


Table 3. Effect'ofi different layer thicknesses. 


v(cm 1 ) 

Nge 

# Layers 

A(v) 

2090.796 

3 

3 

1.12296698 

2090.796 

9 

7 

1.20571646 

2090.796 

35 

18 

1.22967891 

2090.796 

- 

139 

1.231378 

2090. 806 

3 

3 

0.204556963 

2090. 806 

9 

7 

0.202129669 

2090.806 

35 

18 

0.20184368 

2090.806 

- 

139 

0 . 201829374 

2090.816 

3 

3 

0.0613707937 

2090. 816 

9 

7 

0.0616293664 

2090.816 

35 

18 

0.0614016373 

2090.816 

- 

139 

0.0613723083 
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produce thicker layers in the upper portions of the atmosphere where the 
Ui values are small. Since the line center absorbance is more sensitive 
to these upper layers, the loss of detail has more effect on it than wing 
absorbance values. 


The program in Appendix A includes a subroutine called PATHST which 
performs the layer combination. The parameter Nge is used to set the 
fineness of the layers used to calculate the spectrum. As written, how- 
ever, the subroutine combines the layers in such a way that each has an 
airmass at least as great as U/Nge but as close to U/Nge as possible. 

For the spectra in the next section, a value of Nge = 15 was found to be 
acceptable, but Table 3 indicates that different values may be required 
for different types of spectra. 


The convolution portion of the program was tested on a flat 
(ch:/c)v = 0) spectral region with a Gaussian slit function of the .form 


i , N v&r 2 

a Gauss( v - v ) = ^^, exp 


(v * - v)^ fe 2 
H 2 ‘ 


( 26 ) 


A five-point initial trapezoidal evaluation was used to start ttie Romberg 
iterations with a frequency extent of 8H. The T and M table entries cal- 
culated for one run are shown in Table 4. A convergence criterion cor- 
responding to 1$ between the respective T and M table elements [see Eq. 
(22)] was used. Only the diagonal of the T- table is shown since the 
other elements are not calculated in the modified routine. The table ele- 
ments must be multiplied by 2 to get the corresponding transmittance val- 
ues. In this example, the routine went two iterations beyond the initial 
trapezoidal and midpoint evaluations before converging. The final value 
shows excellent agreement with the true value. 


The routine was also tried on a Voigt line with parameters 
v 0 = 2875.90 cm -1 
S 0 = 0.0668 atm -1 cm -1 /cm 
= 0.023176 cm -1 
E" = 400.00 cm -1 
m = 36 .amu 


at a uniform pressure of 1 atm and temperature of 300 K. The total airmass 
was u = 200 atm cm. Table 5 shows the results for 10 different observed 
frequency positions spaced 0.2 cm -1 apart. The second column gives the 
calculated trans m ittance values with e = 0.01, H = 0.083819 cm -1 , and an 
initial five-point trapezoidal evaluation that extended 2-l/2 H on either 
side of the central position. By comparing values equally spaced around 
2875.90 cm -1 the consistency is seen to be ~0.001, which is e/10. The 
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Table 4. A sample T and M table for convolution rou- 
tine. (Table values must be multiplied by 
2 to give transmittance value.) 


TRAP 

MID 

•Iteration 

1 

2 

0.5284401 

0.471553 


•O'. 4999984 

0.4999976 

0.4905153 



0.5094798 

0.4999971 

0.5006297 




0.4993649 

T = 

0.9999946 

(exact = 1.0000000) 


Table 5. Convolution results for a sample Voigt line convolved 
with a Gaussian slit function. 


v (cm 1 ) 

T con 

Iterations 

Total 

Points 

New 

Points 

2875.00 

0.9984500 

1 

9 

9 

2875.'20 

0.9952382 

1 

9 

3 

2875.40 

0.9944796 

1 

9 

4 

2875.60 

0.9885907 

1 

9 

4 

2875.80 

0.85U99 9 

3 

33 

29 

2876.00 

0.8508019 

2 

17 

8 

2876.20 

0.9872358 

1 

9 

3 

2876.40 

0.9942334 

1 

9 

4 

2876.60 

0.9960335 

1 

9 

4 

2876.80 

0.9964244 

1 

9 

4 

TOTAL 72 
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third column lists the number of iterations required for convergence of 
the convolution integral. A value of 1 means that the convergence crite- 
rion was satisfied on the first try. The fourth column lists the number 
of monochromatic transmittances required for the convolution, and the 
fifth column lists the number of new monochromatic transmittance yalues 
calculated at each observed frequency. The increased density of points 
around line center, where the monochromatic transmittance changes most 
rapidly, is apparent. 

Table 6 shows si mi lar results for the identical spectrum except with 
a smaller frequency spacing. The fifth column shows the larger savings 
in calculations for more closely spaced observed frequencies. For refer- 
ence, the 10 points in Table 6 took 0.09 second of cpu time to execute. 

This corresponds to 1.7 ms per monochromatic transmittance calculation. 

A later sample spectrum including 55 atmospheric layers and 13 lines, 
averaged 0.15 ms per monochromatic transmittance calculation on the OSU 
Amdahl Computer. 

Sample spectra for other sets of line parameters produced s imi lar 
results. ;It was found that- extremely narrow and strong lines required 
considerably more iterations than broad lines. This can lead to unrea- 
sonably high numbers of iterations unless some maximum number of iterations 
is specified; usually a value of 7 was used. When convergence does not 
occur on the seventh iteration, the program outputs the most recent ele- 
ment of the T- table, sets an error flag, and proceeds to the next fre- 
quency position. 

A variety of other tests were performed on the program, including the 
effect of different convergence criterion values. More is said about this 
in the next chapter. Overall the behavior was found to be satisfactory 
but one persistent problem did arise; it was found that the convolution ■ 
routine would often converge prematurely. An example is shown in Table 7. 
The rapid decrease in the number of iterations, in a spectral region where 
the transmittance is not changing much, indicates a premature convergence. 
This is verified in the T- and M-tables for these two points shown in 
Table 8. The boxed portions correspond to the values in Table 7 where the 
convergence criterion was 0.005 (or 0.5$). The larger values in Table 8 
correspond to a convergence criterion of 0.0001 (or 0.01$). The premature 
convergence at I 9 OO.O 5 cm -1 is seen to be due to a misleading resemblence 
between the first T- and M- values. Subsequent iterations show the true 
value to be considerably lower than these, which the more stringent con- 
vergence criterion has picked up. This type of behavior is often encoun- 
tered in iterative procedures. This problem was only detected when the 
slit function width was much larger than the line widths, though in prin- 
ciple it could occur in other situations. It was found that the problem 
could usually be avoided by adjusting the convergence criterion and the 
number of points in the first trapezoidal evaluation. The parameters can 
only be set based on experience and knowledge of the spectral region, and 
those lacking both these should indicate some caution. 
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Table 6 . Convolution results for sample line with smaller 
frequency spacing. 


v ( cm 1 ^ 

T con 

Iterations 

Total 

Points 

New 

Points 

2875.70 

0.9694133 

1 

9 

9 

2875.74 

0.9396371 

1 

9 

0 

2875.78 

0.8870856 

2 

17 

10 

2875.82 

0.8118306 

3 

33 

18 

2875.86 

0.7381039 

3 

33 

3 

2875.90 

0.7074896 

3 

33 

3 

2875.94 

0.7380002 

3 

33 

3 

2875.98 

0.8117288 

3 

33 

3 

2876.02 

0.8863926 

3 

33 

3 

2876.06 

0.9391911 

2 

17 

1 

TOTAL 53 


Table 7. An example of a premature convergence. 


v ( cm 1 ) 

T con 

It'erations 

1900.00 

0.96867 

5 

1900.05 

0.99230 

1 
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Table 8. The T- and M- tables for convolutions in Table 7. (Boxed portions correspond to a 
convergence criterion of 0.005. Entire tables for convergence criterion of 0.0001. 
Multiply table values by 0.0625 to get transmittance.) 
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15.87587 15.70601 15.88535 15.83293 15.83359 
15.78528 15.89730 15.82943 15.83364 
15.84218 15.82835 15.83370 
15.83524 , 15.83372 





As a final check on the entire calculation routine, two sample spec- 
tra were calculated under conditions similar to those used by Kyle and 
Goldman. 25 The observer was at a height of 4 km with the sum overhead. 

The spectra for the 1900-1902 cm -1 and 1930-1935 cm -1 regions are shown 
in Figs. 4 and 5j respectively. There are some differences in the condi- 
tions between these spectra and the corresponding high-resolution spectra 
of Kyle. Both have an instrument resolution of 0.1 cm -1 but Kyle used a 
triangular slit function and we used a Gaussian. Also, the values of the 
CO 2 and H 2 O mixing ratios are somewhat different. The calculation routine 
is currently set up for only constant mixing ratios, whereas Kyle used a 
''decreasing Hadmixing ratio with height. An attempt was made to fix the 
COaiHaO ratio in the lowest level. With a CQ 2 mixing ratio of 322 ppm 
(volume) "$ r . this_gave_a^. q ■= JL080 ppm_ ( volume Also,J it'/y/as Ifopuad^that ' 

SNIDER- Will not workjcLth- a; zenith angle., of 0 ° . ^In ~spite-:Of these dif- : 
ferences -the. spectra. show "good, agreement, -indicating the. -two routines . . 
are, 1 consistent . 

t 

In the next chapter, the spectrum calculation routine is used in a 
brief study that illustrates some information analysis techniques for 
experimental design. 



COMPARISON UITH KYLE SPECTPUfl 



spectrum 
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Fig, 5 - Comparison with Kyle 5 spectrum 






III. SOME INFORMATION ANALYSIS TECHNIQUES 


A. IDENTIFYING "BAD" PARAMETERS 

The equations in the previous chapter describe how the appearance of 
an atmospheric absorption spectrum depends on the physical characteristics 
of the atmosphere and the molecular properties of the absorbing gases. 
Assuming that the molecular parameters (line intensity, broadening coef- 
ficient, line position, lower state term value), are known, the equations 
predict what the spectrum will look like, given the atmospheric conditions. 
This report is concerned with the inverse problem (i.e., inferring the 
atmospheric conditions given the spectrum) . It is clear that those as- 
pects of the atmosphere which have no effect on the spectrum can not be 
determined from the spectrum. Similarly, those aspects which have only 
a slight influence on the spectrum will be harder to determine than those 
aspects with a large influence. For example, consider the calculated 
spectrum in Fig. 6. This 200-point spectrum corresponds to the same 
observer-source geometry as in Chapter II-F with the 190-layer model and 
Nge equal to 15 in Eq. (25). The CO2 mixing ratio was 322 ppm, and the 
HaO mixing ratio was 3.889 ppm. Both were assumed independent of height. 
The slit function was Gaussian, with a full width at half height of 0.02 
cm -1 . This particular region contains contributions from nearly 150 CO2 
lines and about a dozen H2O lines, so it should be sensitive to the mixing 
ratios of CO2 and H2O. It contains no information, however, about the 0 3 
mixing ratio. A convergence criterion of 0.5$ was used with an extent of' 
0.05 cm -1 in the convolution. A more stringent criterion of 0.1$ produced 
an rms change in the transmittance values of 2.89 X 10~ 4 with a maximum 
of 0.003, indicating that the accuracy of the convolution is consistent 
with the convergence criterion. A frequency spacing of 2 points per full 
width at half height of the slit function was used. 

Several atmospheric parameters were selected for investigation. 

These are listed in Table 9 along with the two different values used to 
study the influence of the parameter on the spectrum. The spectrum in 
Fig. 6 corresponds to the upper value for each parameter and is referred 
to as the reference spectrum. Eight new spectra were then calculated by 
changing the values of the parameters one at a time. Figure 7 shows the. 
effect on the spectrum of increasing the zenith angle from 95.2° to 95.4°. 
The difference spectrum is shown in Fig. 8. The sum of the squares of the 
differences for this pair is 0.44482, which is listed in column 3 of Table 
9. The corresponding values are also shown for the other seven parameters. 
These numbers are one indicator of how much influence a particular param- 
eter has on the spectrum. Of course larger parameter changes produce 
larger squared differences so the only valid comparisons between parameters 
are conditioned on specific parameter changes; e.g., a 0.2° change in 
zenith angle is more influential than a 2 K change in the temperature at 
layer 4l (18.38 - 18.80 km). It is reasonable, therefore, to scale the 
difference spectrum by dividing it by the size of the parameter change or 
step before squaring. Let 9 i denote the ith parameter with step £ 9 ±. The 
procedure above produces a set of 8 quantities Wn, i = 1, . . . 8 
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Table 9. Parameter values used to compute sample spectra. 


Parameter 

Values 

E (^) 2 

1. CO 2 mixing ratio 

322 ppm 
354.2 ppm 

0.011773 

2. Resolution 

0.02 cm -1 
0.028 cm -1 

0.062913 

3. Zenith angle 

95.2° 

95.4° 

0.44482 

4. Temperature in layer 4l 

216.65 K 

218.65 K 

8.4118 x 10" 5 

5. Temperature in layer 42 

216.65 K 

218.65 K 

5.4798 x 10“ 5 

6. Temperature in layer 40 

216.65 K 

218.65 K 

4.888 x io“ 7 

7. Ground- level pressure 

1013 mb 
1063 mb 

0.008109 

8 . Background value 

1.0 

1.025 

0.10954 
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Fig, 7 - Spectrum for increased zenith angle 95. ^°s 95.2° (reference spectrum) 
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Fig. 8.' (Reference spectrum) increased zenith angle spectrum) 



where 


Wii 




(27) 


with N the number of points in the spectrum, y(^) 

the reference spectrum, and y^ = ^yj^ , . . . , 
the ith parameter stepped by AS A . 


■ fe 0) . 


y, 


i) ) 


N 


the spectrum with 


These scaled and summed square difference values alone, however, are 
insufficient to indicate how well groups of parameters can be determined 
simultaneously. This is illustrated by Fig. 9 . It shows the difference 
spectra for changes in the temperature values T(4l) and T(42) for layers 
4l and 42. From the similarity of these spectra and the assumption that 
the change in the spectrum Zty is linear in /^, it is clear that a simul- 
taneous increase of 2 K in T(4l) and a 2 K decrease in T(42) will produce 
a much smaller difference spectrum than either one in Fig. 9. Each change 
tends to cancel out the other. This corresponds to a 2 K change in 
[T(4l) - T(42)]/2, with the average temperature [T(4l) + T(42)]/2 held 
constant. Similarly, a simultaneous 2 K increase in T(4l) and T(42) will 
produce a much larger difference spectrum than either in Fig. 9. This 
corresponds to a 2 K change in the average temperature, with the difference 
temperature held constant, and indicates that changes in different linear 
combinations of the parameters can have significantly different amounts 
of influence on the spectrum. That, in turn, implies that some combina- 
tions are harder to determine fiom a spectrum than others. Similar ob- 
servations apply to linear combinations of all 8 different parameters. 

The question arises as to whether there exists some linear combina- 
tion 


8 

* : = Z 

j=i 


which has the least influence on the spectrum and hence is the hardest to 
determine. Lees 26 addressed this question in connection with attempts to 
analyze methyl alcohol spectra. He found that changes in a particular 
linear combination of the unknown parameters produced no change in the 
spectrum, which he called a linear dependence, and he proposed a technique 
for identifying these as well as "near" linear dependences. In order to 
adapt this technique to the present problem, form the quantities 


Wik = 




y ( k ))/ae 1 ^ 


(28) 
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and place these in a matrix W, where the ith diagonal element is Wii from 
Eq. (27) and the i - kth and k - ith element is Wq k . This matrix is often 
referred to as the information matrix 27 or the curvature matrix. 28 

To be more precise, W becomes the information matrix in the limit as 

A0 0 and the quantities y^ - y^^/A0^ become the partial derivative^ 
of y with respect to 0 ; evaluated at Go. W is therefore a function of 0 O . 

- ■ - AM 

For small simultaneous changes in the parameters = 1 • I the summed s 

\ a ej 

squared difference spectrum is given by 

Q(A0) = £ £ [4 0) " y k .^)l (29) 

ij k L J 

or, in matrix notation, 

q(3) = 'A9'W& , 

where ' denotes transpose. Lees proposes examining the eigenvalues and 
eigenvectors of a modified information matrix W* given by 

w ik = ^ik/ V^ii^kk • 

Let the -eigenvector A* with elements A* correspond to the smallest eigen- 
value of W*. The corresponding linear combination of parameters is given 

t>y 

8 

0 = £ A? • 9 i • VwTT . (30), 

j=i 

When this eigenvalue is much smaller than the others, he says that it 
indicates a "near" or exact linear dependence. Table 10 lists the infor- 
mation matrix W, and modified matrix W* for the set of 9 spectra already 
introduced. The eigenvalues and eigenvectors of W* were calculated in 
SPEAKEASY with the command EIGENVALS. These are- listed in Table 11. None 
of the eigenvalues is significantly lower than the rest, indicating that 
there are no exact linear dependences among the 8 parameters in Table 9. 

It is still possible to deduce the ,l worst" linear combination from these 
results with the following restriction: 

For each of the 8 parameters 0 £ define a parameter step ££ ± 
such that Q(ds0i) = 1. Assumi ng th at Eq. (29) holds, these 
steps are given by = 1/ Vw 77. The appropriate steps are 
listed in Table 12. Define new parameters {\|r} according to 

ti = 0i/A0i . 
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Table 10. Information matrices 


Mi 


ROW 1 

1.1355E-5 

0.034772 

0.010634 

1.2303E-5 

1.0215E-5 


1.0832E-6 _ 

5.7883E-6 

-0.020427 



ROW 2 

0.034772 3932 

-1.6132 

•0.090516 

-0.072398 


-1.8376E-4 

-0.0091847 

-70.879 



ROW 3 

0.010634 

-1.6132 

11.12 

0.012685 

o.oio46i 


o.ooii462 

0.0059693 

-22.043 



ROW 4 

1.2303E-5 

-0.090516 

0.012685 

2.1029E-5 

1.6957E-5 


1.3054E-6 

7.0995E-6 

-0.021157 



ROW 5 

1.0215E-5 

-0.072398 

0.010461 

1.6957E-5 

1.3699E-5 


1.0741E-6 

5.8427E-6 

-0.017702 



ROW 6 

1.0832E-6 

-1.8376E-4 

o.ooii462 

1.3054E-6 

1.0741E-6 


1.222E-7 

6.2136E-7 

-0.0020892 



ROW 7 

5.7883E-6 

-0.0091847 

0.0059693 

7.0995E-6 

5.8427E-6 


6.2136E-7 

3.2436E-6 

-0.011075 



ROW 8 

-0.020427 

-70.879 

-22.043 

■0.021157 

-0.017702 


-0.0020892 

-0.011075 

175.27 




W*:: 





' 1 

-0.165 0.946 

0.796 

0.819 0.92c 

> 0.954 

-0.458 


l -0.0077 

-0.315 

-0.312 -0.0084 -0.081 

-0.085 


1 - 

0.829 

0.848 0.982 

; 0.994 

-0.499 



l 

0.999 0.814 

0.860 

-0.348 




1 0.83c 

) 0.876 

-0.361 

[SYMMETRIC] 


1 

0.987 

-0.451 





1 

-0.465 


1 
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Table 11. Eigenvalues and associated eigenvectors of y* from Table 10. 


Eigenvectors of W* (as columns): 


ROW 1 

-0.02216 

-0.043359 

0.069408 

0.015821 

-0.25559 

-0.39562 

0.74657 

0.46191 

ROW 2 

0.0038524 

-0.39618 

-0.032105 
0.83294 . . 

-0.052676 

0.06913 

0.18343 

-0.32712 

ROW 3 

-0.008285 

-0.11331 

0.48954 

0.14197 

0.71662 

-0.40684 

-0.13348 

0.17436 

ROW 4 

-0.67639 

-0.048948 

0.08705 

-O.23103 

-0.093432 
-0 . 38268 

0.087113 

-0.56248 

ROW 5 

0.7281 

-0.045501 

-0.032767 

-O.21962 

-0.0074687 

-0.38829 

0.12842 

-0.50117 

ROW 6 

o . 062616 
-0.17077 

0.23679 

0.13061 

-0.59867 

-0.40051 

-0.57696 

0.20514 

ROW 7 

-0.088656 

-0.1211 

-0.83043 

0.068128 

0.22346 

-0.41116 

-0.18346 

0.17485 

ROW 8 

4.4771E-4 

-0.88319 

0.013133 

-0.40311 

0.034914 

0.21564 

0.023464 

0.095041 


Eigenvalue's of : '.y*:‘ 


(1) 

1.0655 x 10 

(2) 

0.0016232 

(3) 

0.010273 

(4) 

0.064734 

(5) 

O.29085 

(6) 

O.68626 

(7) 

1.1985 

(8) 

5.7477 
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Table 12. 

Parameter steps producing a summed 
squared difference spectrum equal 
to 1, assuming Eq. (29) holds. 

(1) 

296.76 ppm (CO 2 mixing ratio) 

(2) 

0.015948 cm -1 (resolution) 

(3) 

O. 29988 0 (zenith angl e ) 

(4) 

218.07 K AT(4l) 

(5) 

270.18 K AT(42) 

(6) 

2860.6 K at(4o) 

(7) 

555.25 mb AP(0) 

(8) 

0.075535 (background) 
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For these new parameters, a unit step in xjr ± produces a 
Q(Aj/) = 1. Now consider the set of al 1 linear combinations of 
steps in the {\Jr} of the form 


such that 


A<i> = I C2&t2 

iCa£tye> 




(31): 

t 


The set of steps A<t> i? with 


• Ci = 1 ; 

1 

and 

Cj = 0 ; i / j , ; 

belongs to this set of A<t> satisfying Eq. (31 ) , and for each Ad>i, ! 

Q(A<t> i ) = 1 . In a sense, steps' Ai> satisfying Eq. (31) are all the steps : 
of the same "size," in terms of parameter changes. From the theorems on! 
quadratic forms 29 it can be shown that the extremal values of Q(At>) for i 
all At> correspond to AVs whose {cj.}, values are the elements of the eigen- 
values of W*. Consequently the smallest Q(A<i>), and hence smallest change 
in the spectrum, for all parameter changes of the same "size" as those in 
Table 12 results from: 1 1 



■0.02216 

X 

296.76 ppm 

0.0038524 

X 

0.015948 cm' 

■0.008285 

X 

0.29988° 

■0.67639 

X 

218.07 K 

0.7281 

X 

270.18 K 

0.062616 

X 

2860.6 K 

■0.088656 

X 

555.25 mb 

4.4771 x 10 -4 

X 

0.075535 


t 


-6.576 ppm j 
0. 00006144 j cm -1 
-0.0024845°! 
-147.5 K 
196.7 K ’ 

v \ 

179.1 k : a ; 

-49.23 mb | 

I 

^ 3.382 x 10" 5 


The most significant aspect of this set of parameter changes is the tern-! 
perature steps. The temperature in level 4l is ‘decreased by 147.5 K while 
the temperatures in the two adjacent , layers are increased by approximately 
the same amount. This zigzagging of! the temperature profile is a problem 
that often occurs in profile retrievals.- 30 The summed squared difference 
spectrum corresponding to this "worst" set of parameter changes is given, 
by the corresponding eigenvalue from) Table 11 ; i.e. , 1.0655 X 10 -4 . It j 
is seen that this is smaller than most of the values in the third column- 
of Table 9. i 


The same quadratic forms theorem allows the identification of the ] 
set of parameter changes A maY leading to the largest change in the spec 7 
trum. It is given by the eigenvector associated with the largest _eigen- ) 
value of W*. From the eighth column, in Table 12 it is seen that A<t> max ! 
corresponds to simultaneous decreases in CO2 mixing ratio, zenith angle,! 
ground level pressure, and all three 1 temperature values. All of these j 
changes tend to decrease the amount of absorption in the spectrum and hence 
add together to increase the summed squared difference spectrum. The net 
change in the spectrum is more than 50,000 times larger than for j 

I 

B. ESTIMATED PARAMETER UNCERTAINTIES ) 


Another useful property of the information matrix W is related to the 
uncertainties in estimates for the values of the unknown parameters to be 
determined from an observed spectrum y 0 t) S . If y(0) is the calculated 
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spectrum for the set. of parameter values 0, then one estimate for the 
true parameter values corresponding to y^g is the well-known least-squares 
estimate 0. This is given by the value of 0 which minimizes 

N - _ 

^yobs,i - • 

-IfL the .rmsiridise .value (is cr, then the asymptotic variance-covariance' ^ 
matrix for 0 is 



where W is evaluated at 0. 

The matrix in Table 10 can therefore be used to estimate, for a given 
noise level, how precisely the eight parameters in Table 9 could be ob- 
tained from a spectrum similar to that in Fig. ;6'. - Attempts to invert W, , 
however, led to warnings by the inversion routine that the matrix is 
nearly singular and the calculated inverse may not be significant. This 
implies that numerical difficulties would be encountered while trying to 
retrieve all eight parameters. If several of the parameters are con- 
strained to given values, this eliminates the corresponding row and column 
from W<q In general this ' m akes r etr i evalV of . the other parameters easier.'' 
The effect of various sets of constrained parameters were studied by 
eliminating selected rows and columns and then trying to invert the re- 
sulting reduced matrix. Some of the results obtained are listed in Table 
13. It is seen from the column labeled 3> that if only the CO 2 mixing 
ratio, resolution, and zenith angle are determined- from a spectrum like 
Fig. 6 with 1$> noise, the respective uncertainties are 10.5 ppm, 0.00018 
cm -1 , and 0.0105°. If the noise is 2$, the uncertainties would be twice 
as large. Any attempts to determine more than one temperature value 
always led to very large uncertainties, as compared to the cases where 
only one temperature was determined. This is consistent with the conclu- 
sions in the previous section. Trying to determine T(4l) and T(42) is 
the same as trying to determine their difference and average. Yet the \ 
temperature difference is very hard to determine, and this gives rise to| 
increased uncertainties in T(4l) and T(42) . 

Another way to think of this problem is to consider the addition of 1 
a constraint as an injection of information into the analysis. The addi T 
tional information spreads out over the variance-covariance matrix and 
leads to smaller uncertainties in the unconstrained parameters. The ex- 
tent to which the constraint of a- 0 1 improves the retrieval of 9j is deter- 
mined by the correlation coefficient between 9i and 0 1# This is given by 

s i)5 /Vs77s77 . 


It is also seen in Table 13 that attempts to retrieve the zenith angle 
and the ground level pressure P(0) produced inversion problems and in-'- 
creased uncertainties. From Table 11 it is seen that this corresponds to 
the second "worst" parameter combination (column 2 of eigenvectors). Its 
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Table 13. Parameter uncertainties for various sets of parameter constraints, assuming an rms noise of 1 
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•Vindicates inversion of nearly singular matrix 










corresponding eigenvalue, though larger than the smallest, is still 5>000 
times smaller than the largest eigenvalue and therefore corresponds to a 
set of parameter steps having a small effect on the spectrum. 

The above techniques of information analysis explore the quality of 
the information in a proposed set of data. They’ .can be used, therefore, 
to compare different experimental designs by comparisons of the informa- 
tion in the data produced by the different designs. This function of the 
techniques was employed to study the effect of decreased instrument res- 
olution on the information in the spectra analyzed above. 

C. LOW RESOLUTION INFORMATION ANALYSIS 

The calculation steps of the previous two sections were repeated 
with the parameter values in Table l4. The new reference spectrum is 
shown in Fig. 10 . This spectrum contains 40 points from 1900 - 1902 cm -1 
which is 2 points per full width at half height of the slit function. 
Clearly, since there are five times more points in the high-resolution 
spectra, the summed squared difference spectra should be correspondingly 
higher. Any deviation from a five-fold change is an indication of dif- 
ferences in the influence of the parameters on the spectrum. The ratios 
of the summed squared difference spectra are listed in column 4 of Table 
l4. Only the value for the resolution is less than the expected value of 
5. Since both resolution changes were +4%, this implies that the low- 
resolution spectrum is more sensitive to percent changes in instrument 
resolution than the high-resolution spectrum. Similarly, the large value 
( 9 . 23 ) for T(4l) indicates that the low-resolution spectrum is less sensi- 
tive to changes in the temperature in layer 4l and hence should yield con- 
siderably larger uncertainties in this parameter. Table 15 shows some 
typical standard deviations for the low-resolution spectrum. It is seen 
that the uncertainties are larger than for the high-resolution spectrum, 
especially the temperature uncertainties. 

Table 16 shows some of the results of a Lees-type search for the 
'bad" parameters. The lowest eigenvalue of the modified information matrix 
corresponds to the zenith angle and P(0) values and the second lowest cor- 
responds to the zigzagging temperature values. This is a reverse in order 
from the high-resolution case for the two "worst" combinations. Once 
again, however, the best combination is a simultaneous decrease in CO 2 
mixing ratio, zenith angle, ground-level pressure, and all three tempera- 
tures. The variation from worst to best corresponds to a change in summed 
squared difference of better than 100', 000. 
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Table 14. Parameter values used to calculate low-resolution sample spectra. (The fourth column 
gives the ratio of the ^ (Ay) 2 values from Table 9 and those in column 3 for each 
parameter . ) 
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Table 15. Parameter uncertainties for various sets of parameter constraints for low-resolution 
spectra, assuming an rms noise of 1%. 
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Table 16. Some eigenvectors and associated eigenvalues of modified 
information matrix for low-resolution spectra. 


(1) 

Eigenvectors; 

(2) 

(8) 

0.0083607 

0.05296 

-0.3934 

0.025936 

0.0074088 

0.11675 

-0.70477 

-0.13512 

-0.39486 

-0.16379 

0.65338 

-0.38916 

-0.17926 

-0.70975 

-0.39096 

-0.016959 

-0.06559 

-0.38969 

0.66577 

0.20956 

-0.39544 

-0.00 71575 

0.0028461 

0.25133 


Eigenvalues: 


(1) 

5.3967 x 10 

(2) 

1.4169 x 10 

(3) 

0.0062167 

(4) 

0.020332 

(5) 

0.12136 

(6) 

0.45162 

(7) 

1.134 

(8) 

6 . 2663 
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IV. CONCLUSIONS 


This report has described the operation of a computer program for 
calculating atmospheric absorption spectra for slant paths through the 
atmosphere, using as input the AFGL line parameter listing. This program, 
unlike many others, was specifically written to employ an efficient means 
of including the effects of instrument spectral response functions. In 
addition the routine is written so that the precision of the convolved or 
degraded spectral values can be varied to match the noise level of observed 
spectra. The program was tested and compared to calculated spectra from 
another source and found to be working correctly. 

Some information analysis techniques were then applied to a particular 
example of slant path spectra. The analysis included an identification of 
the ,f bad" parameters and estimated uncertainties for possible sets of un- 
known parameters. All that was required was a set of spectra and the cor- 
responding conditions for each spectrum. With the aid of a few simple 
matrix operations the character of the information in the spectra was 
quickly determined. These techniques are simple and versatile and could 
be used on many different types of data, both simulated (as in this study) 
or observed. 

The analysis indicated that trying to retrieve closely spaced temper- 
ature profile values can lead to large uncertainties in the temperature 
difference. Problems were also found with simultaneous retrieval of the 
base level pressure and apparent zenith angle. By analyzing a second set 
of spectra at a lower spectral resolution, the loss of information with 
loss of resolution was investigated. In particular, the loss of tempera- 
ture information was found to be more pronounced than that for the other 
spectral parameters considered. 





APPENDIX A 


THE SPECTRUM CALCULATION PROGRAM 

The basic layout of the program is shown in the flow chart (Fig. A-l) . 
The four basic parts are (l) a main program controlling input and output 
and calculation order; (2) the Snider -Goldman slant path program (subrou- 
tine SNIDER and associated subroutines called within SNIDER) for calcu- 
lating the path of radiation in an inhomogeneous atmosphere; (3) subroutine 
PATHST for converting the output of SNIDER into a series of layers with 
specific ainnass, temperature, and pressure values; and ( 4 ) subroutine 
NSPEC for calculating the spectrum. All of the source statements for 
these parts are listed at the end of this appendix. 

I. MAIN PROGRAM 

There are three basic sets of inputs needed by the main or control- 
ling program; these are listed in Table A-l. The first consists of the 
input parameters for the Snider -Goldman routine; these specify the 
observer-source geometry (ALAT, Z, HOBS) and the atmospheric profile 
(LEVELS, H, T, P) . Since SNIDER produces its own pressure profile by in- 
tegrating the hydrostatic equation and using P(0) and the temperature 
profile, only the first pressure value is used. 

The second set consists of parameters needed by PATHST. These are a 
value for Nge, which determines the smallest ainnass in a layer of [total 
airmass/Nge], and a set of values, R(l), R(2), . . . , R(5), for calculat- 
ing the mixing ratios in the layers. As currently written, only constant 
mixing ratios are allowed, but other types can be easily included by suit- 
ably modifying subroutine RMIX. Also, provision has been made for only 
five absorbing species but others could be included. 

The third set of inputs are values needed to control NSPEC. These 
include the upper and lower frequency/' positions, VU and VL; the number of 
observed frequency values, NVO; and the parameters which control the con- 
volution routine. These latter are discussed more fully in Chapter III 
describing subroutine NSPEC. Also included in this set is the unit number 
for the line parameters. This parameter must be passed to NSPEC which 
controls the input of the line parameter data for the spectrum calculation. 

The output from the main program is in two parts. First is the cal- 
culated transmittance values starting from t(VL) up to t(VU). Second is 
a series of error codes. The subroutine ERROR called in the last state- 
ment of the controlling program dumps the stored error code numbers for 
any errors which may have occurred during execution of previous steps. 

These code numbers are in the order in which the errors occurred, and may 
be interpreted by referring to Table A2. No claims for completeness are 
made for this list of errors. An attempt was made to write self- correcting 
routines into the program for some of the more simple errors, but there . 
are undoubtedly many which have not been corrected for. The first output 




Subroutine SNIDER 
and associated 
subroutines 


Subroutines 
PATHST, RMIX 


Subroutine NSPEC 


Fig. A1 - Overall program logic 
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Table Al. Inputs to main program. 


LEVELS 

ALAT 

Z 

HOBS 

NGE 

LUNIT 

H, T, P 
R 

vl; 

NVO 

MM) 

NO ) 


Number of levels in atmospheric profile (must be less than 

190) 

Latitude of the observer (in degrees) 

Apparent zenith angle of incoming ray (in degrees) 

Observer height (in kilometers) 

Integer for structuring layers for input to spectrum calcu- 
lation routine 

Unit number of file containing line parameters. This must 
be the same number as in a JCL DD card; e.g., if LUNIT = 11 
then there is a DD card of the form //GO.FTUEOOl DD . . . 
identifying the line parameter file. 

Atmospheric profile where H is the height in kilometers, T 
is the temperature in Kelvin, and P is the pressure in mb. 
These are on unit 10 (see above). 

An array of mixing ratio values. As currently written it 
requires five values R(l), . . , R(5) for the volume mixing 
ratio of each species; i.e., (molecules of absorber) /(total 
molecules) . 

Lower frequency position in cm -1 
Upper frequency position in cm -1 
Number of observed frequency values 

Set convolution iteration limits (see write-up for subrou- 
tine NSPEC ) 


EPA 


Convergence criterion. This corresponds to the estimated 
uncertainty in the transmittance values. 


EXTNT 

B 


LINTOP 


Sets the integration limits for the convolution. It is 
given as the number of B values (see below) in the integral. 

Resolution width in cm -1 . As currently written this value 
is half of the full width at half max of a Gaussian slit 
function. 

Sets the maximum number of lines contributing to, any one 
spectral value (must be less than 200) 
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Table A1 (continued) 


TTOP Sets the maximum value of the absorption coefficient. When 

a calculated absorption coefficient exceeds TTOP, the cor- 
responding monochromatic transmittance is set to zero. 

NT ST Integer debug flag used to output intermediate convolution 

results. 

NTST = 0 normal operation 

NTST = 1 the M- and T-table values are printed out on 

unit 13 for each step of each convolution 

NTST = 2 the monochromatic transmittance values and 

corresponding frequency positions are printed 
-out on unit l4 ' T 

The NTST = 2 function has been removed from the version 
listed in this appendix:. 

All input is done with unformated READ statements so that there are no 
specific format requirements for the input cards. 


52 





Table A2. Error codes for subroutine ERROR. 


Code No. 

Description of Error and Action Taken 

1 

LEVELS value is greater than 190. LEVELS is set to xninimum 
of 190 and LEVELS. 

2 

SNIDER could not find index of atmospheric layer containing 
the observer. Execution terminated. 

3 

SNIDER could not find index of atmospheric layer containing 
the source. JHT is set equal to LEVELS, where JET is the 
variable labeling the index of the source layer. This is 
equivalent to assuming the source is at the top of the at- 
mosphere. 

-4 

Zenith angle is incorrect for specified observer- source 
geometry; e.g., Z > 90 but HOBS > H souxce . Execution ter- 
minated. 

5 

SNIDER has looped 10 times trying to find astronomical 
zenith angle without success to 0.001. Execution termi- 
nated ,(see Ref'. (10)]. 

6 

GAMMA = 0 in subroutine 0M. Execution continues; resuits 
may be in error. 

7 

Error in OM, not #6 above. Execution continues; results 
may be in error. 

8 

Error in subroutine 0PT1. Execution continues, results may 
be in error. 

9 

FTOT in subroutine PATHST is less than or equal to 0.001, 
indicating that essentially no airmass is in the path. 
Execution continues; results may be in error. 

10 

NM in subroutine NSPEC was specified too large. NM is set 
to minimum of 7 and. NM. 

11 

Failure to converge in convolution routine of subroutine 
NSPEC for current value of frequency position and specified 
maximum number of iterations. Execution continues; current 
value on T- table main diagonal is used as transmittance 
value . 
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Table A2 (continued.) 


Code Wo. 
12 

13 


Description of Error and Action Taken 


NLINES equals zero in subroutine, NSPEC, indicating that 
there are no active lines at the current frequency posi- 
tion. Execution continues; transmittance value is set 
equal to 1. 

The number of active lines, NLINES, is greater than LINTOP, 
the capacity of array IACTV. Execution continues; enough 
lines are kicked out so that NLINES equals LINTOP. 


14 


VU < VL. Execution terminates. 





value from ERROR is the number of type 11 errors which occurred. Since 
this is a common error and may occur many times, only the number is kept 
and not the sequence. Use of NTST = 1 may be used to check specific fre- 
quency positions. 

II. SUBROUTINE SNIDER 

Questions about the operation of this part of the program can be 
answered by referring to Ref. (10). A few minor changes were made in the 
program listed in that reference in order to run it as a subroutine. 

The most important change involved collecting all the error correc- 
tions into the error code listing (Table A-2) . Also the profile input 
section was modified. 

In addition to these changes, statements were inserted to cause 
SNIDER to output the individual layer parameters (T,P,u) as it traced the 
ray through the atmosphere. This is done by calls to subroutine POUT. 

The layer parameters are stored in the array UED, PED, and TED kept in 
blank COMMON. 

An attempt was made to preserve as much of the flexibility of the 
Snider-Goldman program as possible. This includes the option of specifying 
the astronomical zenith angle and having the program find the corresponding 
apparent zenith angle by an iterative search. To use this facility, the 
input to SNIDER must be changed so that LSTOP, the 4th input value passed 
to SNIDER, is equal to 1; see Ref. (10) before using this option. 

One minor problem encountered with SNIDER during debugging is that 
it will not run with an input zenith angle of zero. The associated sub- 
routines are OM, 0PT1, ASINX, ZCQMP, QUAD, GETT, and POUT. 

III. SUBROUTINE PATHST 

This subroutine takes the individual layer parameters stored in UED, 
PED, and TED and rearranges them into a set of layers with nearly the same 
airmass values stretching from the lowest layer crossed by the ray up to 
the highest level. 

First it sums all the airmass values to get FTOT, the total airmass. 

It then folds the parameter arrays to combine contributions from the same 
layer that occurred during the descending and ascending portions of the 
path, if the zenith angle is greater than 90° . It next starts from the 
lowest layer and combines the contributions of the next highest layers 
until the summed airmass exceeds FTOT/Nge. A call to subroutine RMIX re- 
turns a set of mixing ratio values for this composite layer and then the 
program begins combining the next layers until the airmass again exceeds 
FTOT/Nge. The recombination continues throughout the entire path. The 
resulting modified set of atmospheric path parameters is stored in U, 

PRES, TMP, and RMIXR in blank COMMON. 
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IV. SUBROUTINE NSPEC 


This subroutine takes the path parameters from PATHST, the line pa- 
rameters from unit LUNIT, and the frequency and convolution parameters 
from the SUBROUTINE call and calculates the transmittance. The file of 
line parameters has been created by a separate program BOB listed in 
Appendix B. These consist of (l) the lowest frequency at which the line 
has significant absorption, (2) the line intensity, (3) the pressure- 
broadening coefficient, (4) the lower state term value, (5) the molecule 
code number, and (6) l/2 the frequency range of significant absorption by 
the line. The program has a set of arrays in which the parameter values 
of the active lines are stored. The input file is sorted on the lowest 
contributing frequency values and NSPEC proceeds from the lowest observed 
frequency value, VL, to highest value, VU. The input file is processed 
in sequence and the active lines are sorted by NSPEC according to the 
highest frequency at which they have significant absorption. Because of 
this double sorting (i.e., sorted input file and sorted active list), 
when the transmittance at a, new frequency position is to be calculated, 
the program needs to check only the first active line to see if any lines 
need to be removed from the active file. If the first line is still 
active at the new frequency, all the other lines are also still active. 

This eliminates the need to check every active line at every frequency 
position. If the first line is no longer active it is deactivated and 
the next line is checked. This continues until all lines which are not 
active at the new frequency position have been removed. The program then 
checks the next line in the input file. If it is not active at the new 
frequency, none of the lines below it in the file will be, and the program 
proceeds to calculate the transmittance. If the new line is active, it 
is inserted into the active list in the proper sequence and then- the next 
input line is checked. This minimizes the number of checking operations 
required to maintain a correct list of active lines. Also, by not storing 
the parameters of inactive lines in core, the overall storage requirements 
of the program are reduced considerably. 

In addition to the input parameters, the program stores the molecular 
masses (in amu) in the array FMAS. This value is needed to calculate the 
Doppler width. The active line parameters are stored in the arrays VO, 

SO, AMAS, ALL, EPP, MOLNM, and CENTW. Instead of sorting the elements of 
these seven arrays, an array of indices is kept, called IACTV, and these 
are kept in such a way that the order of the elements in IACTV determines 
the order in which the lines have to be checked for possible deactivation. 

In order to calculate the convolved transmittance, T con , at observed 
frequency position, v*, the program must calculate the monochromatic 
transmittance, T mon , at frequency positions vf, . . . , v„ (refer to Fig. 
2). These are then- multiplied by an instrument response (or ’slit-) func- - 
tion value and summed to give T con . The storage of the T mon values in 
array T proved to be the most logically involved aspect of the calcula- 
tion. The actual values are calculated quite simply using the Hui et al. 23 
Voigt routine and summing absorption coefficients by layers and then by 
active lines in two nested DO-loops. The slit function value, called WT, 
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is evaluated, multiplied by the T mon values, and summed in another loop. 

The problem is how to store the T mon values so as to be able to retrieve 
them efficiently. This is al 1 done by proper arrangement of the values 
in the array T. 

First, the maximum number of T mon values is determined from NO and 
NM according to NTOT = 2 • 2^ + 1. The % on values will then be stored 
in the first NTOT elanents of T as‘ : if each position in T were a frequency 
position with a spacing between elements of DELVC = EXTNT * H/ 2^ + -*- . This 
spacing fits the NTOT values into a total frequency range of EXTNT X H, 
where H is one-half the full width at half height of the slit function. 

To start off the calculation, the first element of the T-table is 
evaluated using NO to determine the number of points, as shown in Table 
A— 3 . Also listed are the number of new points (i.e., T mon values) used 
in subsequent iterations. The index n along the left-hand side of the 
table is used by the program to determine the maximum number of iterations. 
When the n value exceeds NM, the convolution stops iterating and takes 
the convergence failure error route. 

For example, when NO = 1 and NM = 5 > the first trapezoidal rule 
evaluation is made with five points. The corresponding x mon values axe 
loaded into T. The next step involves a four-point midpoint rule evalua- 
tion. These two values are checked for convergence. If the convergence 
criterion is satisfied (AMTS.LE.EPA) , the convolved transmittance value 
is loaded into the array TCI. If there is no convergence, an eight-point 
midpoint rule evaluation is done and the T- and M-table values calculated 
and compared. The procedure continues until convergence or until the 
32-point midpoint rule value is calculated. After this value, the con- 
vergence failure route is taken. In this way, up to 65 T mon values would 
be calculated and loaded into array T in order of increasing frequency 
position. When moving to a new observed frequency position (next VI 
value) , the new i: mon values that need to be calculated are written on top 
of the old x mon values which are no longer needed. In this way space is 
conserved and the values of T can be efficiently retrieved. 

As currently written, NSPEC generates its own observed frequency 
positions at which to calculate T con values. This can be easily changed 
to calculate x con values at specified nonuniformly spaced frequencies by 
passing the array of frequency positions to NSPEC and changing the-scalar 
VI to an array V(l) of these frequency positions. • 
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00002890 
-000 02 9 00 - 
00002910 
-0000-2-9-20- 
00002930 


(— J T 

00002950 


L- / W ^ 

00002970 

A A A A A A n A 


’ KJ S V ’ 

00002990 


RcFk =<DELZ +SECT1) 

00003000 
. nnn n i n 

in l ! |\ j ■” K L 1 K ■ K A U j C U 

2010 DELRS=REFRS 

00003020 

zrtTtjrJ i i\ UU J ! 1 i z. 1 — 

IHTi = IhT + 1. 

If! T i'iCt i ! _ i T _ T i 1 T i ) ftP — — L17 -, 

— — “ " ' uoUU‘3 U JU “ 

00003040 

nnnniimn 


1360 CONTINUE 

DU 13 - 7 5 I=IHT - 1tIND 0 

U = FCNU( FC.MA( ESTR ( I ) ) , FCNB ( E STR ( I ) ) , E STR< I ) , QMSTR A ( I ) ) 


00003060 

000 04 0-70- 

00003080 
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DtLZT=DcLZT+DELZ 00C03100 

V A RC0 M -( ' l 9 ) -RHC < I ~ ) OQ00 3irlO- 

CALL w‘U AO ( TrtET AS ( I + i ) ♦ Z ( I) » H ( I ) , T ( I ) , P ( I ) , GAMMA ( I ) , ESTR ( I ) , ADD 1 , OOOC31 20 


CALL P Ou t ( A uD2,ADD3,ACD4 ) 


00003135 


SECT3 = SEC T 3 + ACD3 


00003150 


DECK =(iJcLZ+AL)01) 


00003170 


R£FKS-RtFR*RAUSEC 


1 37 5 CUNT Ii'jUt 


RULAST=kHG(JHT-1 ) 


00003190 
-eoee- 3 - 2 - o &- 
000032 10 
-0GQO3-2-2O- 
000032 30 


ELAST=ESTR ( jHT-1 ) 


000032 50 


Tl A S T= T (JHT-1 ) 


HL A S T-H ( JriT- 1 ) 


00003270 

0000-3-2-80- 

00003290 


CALL UM(HGHT,HLAST»H( JHT-1 ) ,TLAST,GLAST,OLAST, XXX) 000033 10 

U=FCN(j( F C ,\i A ( ElAST) y P C NOT E l - A - ST - ) t E L AS T , OL A S ~ T ~) 0OOO33-2Q- 

CALL aCuMP(ZLaST,ZFIN,HLAST,HGHT,L,DEL7,THFIN) 00003330 

R 0 L A Th T — QOOO-3-3-AO- 

CALL *UAO( THFIN, ZlAST*HLAST,TLAST ,PLAST,GLAST,ELAST» ADD 1 tADD2* ADD300003350 


CALL POUT (ADD2»ADD3»ADD4) 


00003365 


SECT4 = SECT4 + ADD4 


RcFR=DELR+REFR 


00003380 

-00003-390- — 

00003400 


StCT2=SECT2+A0D2 


IF ( INS. EG. 2 ) GO TO 1355 


GC TO ( 39^0, 1355, 1356, 1358 ), INS 


00003420 
-0000-34-30— 
00003440 
-00003 4-50— 
00003460 


ZLAST=ZO 


00003480 


HLAST=OHGHT 


00003500 



00003520 


REFR=u . 


3950 CONTINUE 


00003540 
-00-00-3 5 50- 
00003560 


SECT4 = SECT4 + TSAVE 


00003580 


SECT2=S£CT2+OSAVE 


00003600 


Z£Nu=ZOS+DERF 


TEFF - 5ECT4/5ECT2 


00003620 

-003036-30 

00003640 


3 9t> 0 CONTINUE 


00003660 


3990 CONTINUE 


00003680 


6k 



3991 

FORMAT ( 1H1»5X25HRE FRACTION AIR^ASS TABLE S , 5 X 1 3 A6 , / , 


00C037 00 


j^UXjoHhLL He IbH 1 j AKt I , n G t U ri b 1 R i L n t i c r\ o / ? AUa» 

3ES ARE I N DEGREES KE L V I N • / , 9CX * * A LL PRESSURES ARE IN 

h L L 1 C’-r CttM 1 U^U UUUj 1 L U 

N/SQ M V 00003720 


MHOa^UHThL A"PP AKcN I A 9 u A l R uNU” I G AL ANbLto AKtT7 
54GX34nIW DECIMAL DEGREES FROM THE ZENITH/) 


UUUUJf JU 

000037 A0 

— c 

3993 

— Wts l i c — t — 3 — f — j -j -j-o — ) — w A V t: L 

FORMAT ( 9CX • TriE wavelength for this table IS ' * F 5 • 

2, • 

MICROMETER 000037 60 

c 

i s-. / / I 
WKITE< 6,3992 ) 


UUUU J 1 /tJ 

3 992 

C 

3-995- 

' PO K M A T '( i I a 5 H L u ri t R"t j~a 'JHU P P fcrR * 3 X 4 H *1 IN • * 7 "X’7 HC H A P n A N * 

W k I T E ( 6 *3993 ) 

--Ft) R ,9 A -T-H X6-HAPP. ~Z~r9-X 6 H H£-H? H-T~t 9 X 6 H F E I &HT , 9'X6nHEi:-&HT-,-5 
18HAIR MASS ,9X5HHDIFF , 3X 1 0HR eFR ACT ION* 9X6HAST. Z,6X 

j a mu rm. a lj u u ltu j f to 

X-8-HA-IR NA S S~, -9 X-00003-8 TO 

6HP (EFF ) ,9X6H00003820 

C 

ri~t err) , 6XOHP ( T AiV) 3 

Wk I T E ( 6 , 3 996 ) 


UUOUOOO U 

000038 A0 

jyV'IT 

ru Pv 1 1 M 1 l JAJ n l 1 j f 1 Ajnt ^ } * 1 A JTl l 5 J t | A / * O A S> H \ J * O A DFl K O 1 

1 6X3H{ 3) t9X3H(9) *9X4H( 10) *6X4H( 11 )* 8X4H(12) */,lX ) 

T TOX Jn'rrn u'O 0 ©"3 0513 

00003860 

3999 

tTn C^TC 

LINC = L INC- 1 


UUUUjU (U 

00003880 


tr IN — - reiver^ — KAUCuN 

QHGHT = FZC ( 0HGHT ) 


UUU03 8 90 
00003900 


IIGHT ~ rZC ( HuHT J 

HMNA = FZC (HMINA) 


UUU0'39IU 

00003920 


o L ALnP — = — K"S“r A K * — 1 r / A I Kto'i — / — b 

xp = ( RADIUS + 0HGHT ) / SCALHP 


UUUU393U 
000039 A0 


"Wtrr — = — StrRr — ( — u «5 — * — JOr — ) — *~ -A B S — (COS ('Z'F-f-Ni - ) 
FERF = EXP ( (WHY ) *-2.0 ) - ( 1»C - ERF (WHY)) 


00003950 ‘ 
00003960 

“TUH O 

9150 

L u N i 1 w u r — 

CONTINUE 


UUUUH ZiU 

D000A270 


L J 1 U r J 

GO TO (1250,9100,1250,1256,1000) ,JMP 


UUUU f ZJ 

00009280 

i luU 

Curt li iM u n 

IF ( ABS(ZT-EEND) .LT. 0.0001 ) GOTO 9200 


■ U 000 ^'2' 90 


tU J“ c 1 u L U 1 

Z0 = Z03 

Cf. TO 1 9 lj L ..... _ 


UUUU H’J”rO 

00004320 

9200 

J M P = 1 

GO TO 3930 


00004340 
n p 004 3 SO 

4 500 
a s 1 0 

CALL ERROR ( 9 , 0 ) 

p. . o v* a t / i </ r i /, c _ Arv-nu *v«v>v \*n cm i nriM \ 


ODDP/t? 7H 


LINC = l INC- 1 

P, CTl | O M 


00009380 

9600 

RtFR S = C . 
SECT ~>-C 


00004400 

00004410 


S6CT3=0. 

— cpf-T/.-n 


00004420 
nnonv, ^ o 

. is. 7 n n. 

GO TO 3950 

lie ( a »-•. c_ f _7__t. i c \,n » i t n n i i r n r n /, on n 


00004440 

000 04 U SO 


L I NC = l I .mC - 2 


00009960 

9 705 

/ n n "3 

FoRo AT ( 1 9 ,37h***i\0 GC0L SOLUTION FOR LINE ABOVE***) 

c . . O A., A t / i * * ,iMOrro*rTrn .dav c t o t •/ c c n t c r 


000095 CO 

• \ rs r\ rs r\ /. c: -> r\ 

1 A C 0 

RETURN 


^ 

00004 S 4 n 


RET ORi'i 

c i n 


noon/. 


S JBkOuT INE U"l ( H,HG ,HNI,TN, GAMMA ,GM EGA, TO ) 

- £.r. .-4iu r> f A-&1T / r* _ a t n ijt C t a d o a n r i C r hmc t i 


00C046 10 

0 0 0 O Z, P ^0 


GM=ABS ( GAMMA ) 

TP ( . -A } T P r, O r : nr. n ■» i p rp HOI 


00004630 

O O O O 4 A 4 O 


0NEb= i •G-uAMMA 

— Gw o ~ Qri EG — — — “ — 


00009650 

0G0 096-60 
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IF (RAC IuS+Kn 5 500 t 500 »t> 0C 

00004670 

*3 00 

0 U "J 1 jTN 0 c ~ — 

X = CGMSTl-h-( RAD I US -HC ) /TN 

OOOUHOoO 

00004690 

600 

o u T u 7 uC 
CONTINUE 

0G004700 
000047 10 


“A - L u N :> 1 J. v ( ri — rrtn 

T0 = TN + CCNSTi*( HO-nN) * ( 0NG) /GAMMA 

UUUU4/2U 
000047 30 

7 0C 

X — A / i u 
CONTINUE 

00C04740 

00004750 


I -* — A / JMIV'IA 

IF ( GNtG.E«i.u. ) GO TO 600 

00004760 

00004770 

4u0 

i i l JL • U Unco* Y ) OUjfoJjf tUU 

I F ( l.o- ( 1.0-0NEG*Y ) ) 8 0G » 800 « 300 

UUU U^ i oU 

000047 90 

JUU 

Uii c bA-ii tTi unTvrnr 1 v \ 1 • /U'NLb ) 

return 

“UUU UHo uu 

000048 10 

ouu 

1 

l iiyuL 

IF( Y.lT. (-30.0) ) Y= 0 . 0 

cno«i at ( i i-| c;o* p\ 

uUUU'to <l'U 
00004830 

1 

r u rv n A i 1 irfi f Dec tnrtrr 

I F ( ABS(Y ) .GT. 25.0 ) GO TO 1C 

UUUuHU MU 

000048 50 

10 

--Url'trb A — i • U Cat l Y“J 

CONTINUE 
.ncT iPi* 

UUUUhc cU 

00004870 

801 
Q t . p 

C AL^ c R kOR ( 6 * C ) 

UUUU4U uU 

803 

r U K F 1 A 1 — \ L A v On — OAflnA-t 1 1 

call E RR0R ( 7 i 0 ) 
p . i n a t / • • > nr, p n n * /, r on o » 

UUUU4VLU 


ru/\ rr A 1 \ dUUnuo f *tt rJ * o J 

RETURN 

C .Vi r\ . . _ 

— — — trU Ul)TTZ'u 

00004950 
nnonzi-O f- n 


subroutine opti ( ho * h *to» gamma* rhqc:* omega* a irmj 

UUUUTvOU 

00004970 


“tro rTTTtjri / jL A pj t~/ b f A'rK'VT 1 1 AKf KA Uru t—± 

RA0-RAD IUb-hO 

" UUUUMVOU 

00004990 


“~a1\ -L,UiM j'I 1 v K'M'U'r l 1 U b A'n'n A ) 

SmALLH=(H-HC)/< RADIUS- HO) 

UUUUjUUU 

000050 10 


u n l i i — l • u o i i allii ~ 

On EG = 0 A MM A - 1 .0 

UU JU lU 

00005030 
nnnn ah lc\ 


”"Af\ c - 1 i *\J — "U i n t b o) II A'L tn r A l\ # / U - 4 lT^ 

S u ‘M U - 0 • o 

Cl I -* V/ -p _n ... . . 

U UU U DUHU 

00005050 
nnnrmoA p 


o Url v — U • U 

XN = 0 .0 

00005070 
onnn^pi q n 


u l h u n - l • j ■■ 

TERMU= 1 .C 

UUUUDU UU 

00005090 
- nnnr^ inn 


It. \\ I'l V — I • u 

TEST =.0 GOG 0j01/(XK*GAMMA ) 

y r,| j T - O . 

00005110 
nnnn r i 

1000 

CONTINUE 

00005130 


m *\ UU'N 1 * 1 11 

X N = X N + 1 .0 

UUUUD1 HU 

00005150 


t7t_ .MU 1 — JLiiU |TUl\Lb 

FmCT=aN/ (XK*DENGM ) 

UUUU^l OU 

000051 70 


1 LkMJ- 1 tKAU'-'l Al 1 

TERM V= TERM V*FACT*XK2 

"UUUUD1 OU 

00005190 


tro 1 u = j U i*! J * i n 1 , n u — 

SUMY=SUMV+TERMV 

UUUU^6 uu 
000052 10 

1040 

ri \ t\UUN 1 zriri i O 0 j » 1 O o 0 

IF ( m 3 3 ( TEST) -A3S l TERM V ) ) 1 0C0 , 1 C0C , 1 050 

— - — 

0C0052 30 

— 1 u > u 

1075 
1060 

“it l Auyt l‘- j 1 J AD'5 t 1 LK 10 ) * I J UUf L U ‘ J L j i ' J I o 

IF (K0UNT-3G) 1100,1100,1080 

WRITE ( 6-y ICOj ) 

UUUUDl HU 

000052 50 

900-05-2-fce 


1065 

a r> 


FORMAT ( 1 7n0 

a r o _i_n__A_ri r i_< C_- 


80GB00 IN 

>'»*•? / n t r\ \ n l 


CRT 1 ) 




00005270 

O r\ on coon 


J 1 



FUn!CTIO.< a S I N X ( X 


ASIimX = X*( 1 • + Z / 6 « + o 07 5 * Z ^ * 2 ) 


000053 10 

-00-005-3-20 

00005330 


I * U <f I . 

COfViOiN. /SLAr.T/ G,AlR/JT ? RST^R,RACILS f CONSTl 

S^ iAL Ln = ( IIG - HO )/( RADIUS - MO) 

ON £ n = 1 • C— S f-iM L L H 


00005370 


00005390 


C-CuS ( Zu ) 


00005410 


CuNE =C #0 I e H 


DEL 2 =S* l US-TWCH ) 


00005430 

-000-054-40 

000054 50 


Otl=S*TwOH/(CCNt+SQRT( CPLUS ) ) 


ZS=ZO+DcLZ 


00005470 

- 0000 - 54-80 

000054 90 


Rt TuRi>< 


000055 10 


SUbROUT INc uU AD ( ThETA 9 ZC 9 HO , T 9 P » G AMM A , FC » SUM t V ALM , P B AR 9 1 BA R ) 


DIMENSION VARCOM 50 ) 


DIMENSION XPTS(6) ,THTPTS(6) tWGHT( 6) t VALUES(6) 


00005530 

00005550 

- 000 - 6 - 55 - 60 - 

00005570 


EQUIVALENCE ( VARCOMC 1 ) »G) , ( VARCOM ( 2 ) ♦ AIRWT ) 9 (VARC 0 M{ 3>9 

-t RS- T -A- R - n ( V AR - e0M - ( ^r)^-R A - ei t>-Si-9 — ( V AR - CGm5 -) 9 CONS - T - 1 )-Ti-V- ARC0 M i-6-)- 

2 9 XPTS (1 ) ) , (VARCOMi 12 ) 9 WGHT ( 1 ) ) 9 (VARCOM ( 18 ) 9 FAKE ) 


00005590 
-00005 600 - — 


EQUIVALENCE ( VARCOM* 19 ) 9 RHO ! 


- 0000 - 5 - 6 - 20 - 

00005630 




■V A*JLVC V 1 PJK •WAV ■ 


l)/( 1 .U+UOGOFL ) **2 00005660 

F - CNu (u OOirL) = 2 .0*G00irL* ( 1 .0 * 0001 TL/ 2 -^ 0 * 000056-70- — 

FCNU ( G002FL, J003FL 9 QOO 4 FL 9 000 5FL ) = 2 o 0*Q004F L*Q00 5FL*Q002 FL / < 1 . 0 +Q00000 5 6 80 

. 03 F L -»~Q0u 5F L -j 000056-90 

FCNh (wOG3FL, Q009FL ) 2 . 0*S I N ( QOO 8F L/ 2 . C ) ** 2 + COT AN < Q009FL ) *S IN ( Q008F000057 00 


SmaLLR=kAO I uS+HO+hO** 2 /(RADIUS-HO ) 


00005720 


A=FCNA ( EO ) 


Du luQO 1 = 196 


000057^0 

-00005-7-50- 

00005760 


HINT (I )=FCfJn(THTPTS( I ) 9 ZO ) 


00005780 


DU 2 00 G 1 = 196 


00005800 


* _ n I 1 M M V ) 
j UUrt n r 


U=FCNu ( A 9 b 9 EO 9 OMEGA 


000058 20 


VALUES ( I ) = 1 ./ SORT ( ARG 1 ) 


00005840 


2000 CUNT Ii'iUt 


00005860 


VALh = G • 


C 0005880 


SUM= VaLuEj ( I ) *wohT ( I ) + S U M 


3000 CONTINUE 


00005900 

-QTTQ-Q5-9-1-0 

00005920 
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V ALf.=v AlCSVmLK 

000059^0 


r u A k — V AL.T’-r 
Td Ak — V A i_ 'V-' T 

“ DutrUJTDt) 

o 0005960 


R l T uRi^ 

E,\iD 

00003970 

00005980 


buLirndU I rYt ut i i ( Hi tn^ f r i » ! i t 1 / j - """ 

CuMNO.M /SlAuT/ OtAlRWT tRSTARf R AD I LS * CONS T I 

00005990 

00006CC0 


1 1“ ( 1 2 • b i • U • J u u rU I f U tj 

K = C 

UUUUou 1U 
00006020 


X — r c / r 1 

XL=mLOG( X) 

000 OfrO‘3 0 
000060 AO 


Ai-L JuoTl''- Ino” HI ) / i 1 
CRIT=tXP(-Al ) 

000060 50 
00006060 


Al=i»/AL ” - ■ ■ 

IF ( AOi ( X-cRIT ) .uT. .OOGO l 1 ) GG TO 100C 

0 0 0 060 TO 

00006C 80 


T 2 - T i 
RcTuRN 

00006090 
0 0006100 

1 o uu 

a C - A~i vXL^mLuG l A i ) 
ACR 1 T= - 1 • / XL 

00006 no 

000061 20 


A2 = • OOOoO i 1 

IP(Al.LT.ACRIT) A2=A2+ACRIT 

000061 30 
000061 40 

1 1 00" 

r — AZ^^AL+’ALojo ( A 2 ) ■” X C 
FP = XL+ 1./A2 

000061 50 
0C0061 60 


A 3- h r / r P ' 

I F ( N • LT • 3 ) GO TO 1200 

000061 70 
00006180 

1200 

I P ( A ( AS- A2 ) •LI •• 00 OOOOlJ 00 TO 1300 

A2 = A3 

u 0 0 0 6 190 

000062 00 


A = ?\"+ l 

GU TO 1100 

UUUUu£ LU 
000062 20 

i jUu 

~ ** — — 

RET URN 

UuUtrrrzr ju 

00006240 

i 1 ou 

i t" \ j t • t y '* 1 c. J o U I U i 0 00 

P2 = P 1*( T2/TI ) **( < CONST 1*( H2-HI ) ) / IT1-T2 ) ) 

LUUUOi DU 
00006260 

1800 

k l i — . 

P^ = P 1* EaP ( C0NST1* ( HI -H2 ) /T 1 ) 

ijUUUCJZ! fU 

00006280 
CinnnA o on 


END 

cimuni it i m c oni it i ai . & ? . s ? i . 

00006300 

nnnnA -xnn 

0 1 Mb NS ION U£D(3GC)tTED(303)fPED(3G0)?H(200)jP(2 00)tT(2 0G) 
r ii m ! v| n a m . p . t . it p n . p ^ n t t c n . t *= n 

IF (Al.EQ.-l.C) GO TO 100 

UcD( IdO) =A1*G.215335E26 

p f-n ( Tr n ) -i . ... _ 

TEDl ItD )=A 3*0.215335626 

100 

IcD = C 

■O-w.T-1 IDs, ... ... .. . .. .... .. 


END 

c i zj n , t r vi r c_atuv o\i j c w t n n \ 


Di-'icNS IGH JIN ( 300) »PU ( 300) »TU( 3C0 ) »U( 50 ,PRES( 80) , TMP( 50 ) , 

'1 1. ^ tv ■ / r. - o 1 t_. t / r\ 1 o o f n r\ > t rw o nr\ \ r> m # c \ 


C Li i v it*i 0 in tlT»PPjTP*UiN»PUtTU»IED»U»°RES?T w ip,RMiXR 



DU 30 1 = 1 , IcO 

PTPT- p-T- i T i 1 ! T M / T \ 



j l U 1 1 l u J * Ul Tl 1 J 

OG J.0G I = 2 ♦ I ED 


100 

ii ntrjtrrtnTrnri u rivl ft i • L _ • 1 • J n 1 i — L i > i u — irJ ■ 

C u N T 1 1 n U t 
rni n-n 


bcQ- 

In r UL u - u 

GG TO 220 

-fVEOirB - I - I 
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Du 200 j=l»NFGL0 

K-NFOLO+ J 

UINlK)=LIrt(J)+UIN(K) 

Pu(i\)-?u(J ) r"b ( K ) 

Tu(K)-Tu( o l + Tll(K) 

■2tr0 — CUNT I,;Ut 
2^3 Jl = .\F0Lu+l 

Kr _ Q 

IF IFTOT.LE.O.Oul) CAlL ERRCR(9,0i 

Su'+u-'crra 

Su'-1PU = 0.0 

Su,MTU= 0.0 

Ou h Ou J = J 1 1 1 EC 

SUmU-S Ui v i U * U IN ( J ) 

S U Pi P U = S u M F U + P U ( J ) 

TO t TlH - j ) 

IF ( SUFlU.GT. FTCT/NGE ) GO TO 390 

G G— - T O V O - 9 

390 K = K + 1 

U ( K ) -5Ur -i U 

PRES ( K ) = SutfPU/ SuMU 

T - nfi -F K - )— ■ 

SUMu = 0 .0 

S u M P U = 0 « 0 

SUMTU=0 « 0 

IF ( J.Eu.IED ) GO TO 420 

900 CONTINUE 



U ( K ) -SUMU 

PkES(K ' )-SUH P U/SU M U 

TMP ( K ) = SUi v iTU/ SUMU 

- 9 - 2 - e --e0fTT-m-i3E 

00 oOO J = 1 » K 
Eh3 5 Ou I J - 1 , 5 

RM I X R ( I J , j ) = RMIX(IJ»PRES(J ) , RN ) 

-5- 6 ~ 0~ CG N TT f\ ‘ U v — - 

600 CONTINUE 

LEVTOP = K 

RETURN 

Ett3 


FUNCTION *MIX(I*P,RN) 

— B - I - M EN -S- IuN RISt ( 5 ) 

RMIX=RN( I ) 

R C T uRiv 
END 

SuOfsuuT I ‘IE i m 3 P EC ( V L » VU , N V Q , iF-’-» NO » EPA»EXTNT,H,LINTGP»NLEVLS, 
*TTGP,NTiT, IuNiT,TCI) 

~ ' C b M P L E X * 1 6 A 0 / ( 1 . 2 2 6 0 7 9 3 1 7 7 7 I Q^ r Ert&Q 2 j 0 . -0 00 - ) - / - , - Al - / 

-M2. 1436 2 36 869970692 50 02,0. ODD) / ,A2/( 1 . 8 I 9 2« 53309 2 1 8 1 5490 02 , 

+ QvOO C ) /,Aj/l 9 . 3 15 ir 56 04 561 38441001 ,0. 000) /,A4/ 

+ (:>. ulo0i4219 62 lu58 9L) 01, C. CD C)/,A5/(5. 9126 262 09773153 00, 0.0 DO)/ 

* , A 6/ 1 u » 5 6 — " I o95 ' 8 35b2ol 5 ' 0 TTT& . 0 DC r 6~ / T~ 30 / ■ ( ■ 1 . 2 2 60 7 9 3 1 2 , 

♦0.Gu0)/»ul/{ 3. 52730625110963558002,0.000 )/• 82/ 

+ ( n . 5 7^ 3 h 47 07 - 03 S 977 3 7 0 G 2-y0 - . 0 00)/, 33/(3. 40 7 1 9-4 G 57 9 2 D 02 , 

+0.000 ) /,o4/( 1 .70354001821091^72002,0.000 )/, 85/ 

+ ( ■ 5'.399 - 290u912 94G2C7D01 -r 0.00C) -- /,D6/ -l- T . 04 - 7 - 9057 1 1 4 26 - 03 9- 900 - 1 « 

+ 0.000) / ,Zri»i-,0C/iPLX 

01 fe N-3 IJN HA-f eOO ) , 0 X f 200 H TX ("2 OC ) ,U I NX ( 3 060 - , PUX ( 300 ) , T - U ~ X - ( 300 ) 

REAL-0 ZVuIGT, X , Y, A8SK , VO, Vi , VK, V IK , VL . VU, V8H, VONXT, VOL, VUP,D2A8S 
DIMENSION VO 1 200 ) , SC < 2QJ ) ■ , ALL ( 200 ) , E P P ( 2 0G- >- , T-C I ( 3 00 ) -, 
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500 l*ll=.Mll + l 

RET U R n — 

1000 w k I T E { 6 * * ) i*i 1 1 » I E R 
R r i U Ri n 
Ei\D 

/* — 

//GL.SYSLIN CD GS.\i=TSC462 .NSPEC .OBJ » D ISP = SHR 

//GU.FT1 JF001 trO DSiN -T SC 4 6 2 - . S M C E R . D A T -*r(-S7W-E' ) , 01 SP-SHR - 

//Gu.FTllFOOl u3 DSK=T SG<t62 . L I MP AR « DA T A , 01 SP = SHR 

// Gl . i T - l -2- r - e— 1 c-8 OYt& O' T- A 

//Gu • FT 1 3r 00 1 00 SYSCUT^A 

- / / GO «'^YSIi\ &0 ■* 

190 95.0 95.0 h 5 . 0 15 11 

— 3 « d ~9~E — &■& — 32 2.0c — SYb — 0.0 — O ' . 0 — 0 ~ o 0 

1900. 0 1900. C5 271 0.C05 5.5 0.05 20C 10.0 0 

■/* 

// 
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V 


APPENDIX B 


THE LINE PARAMETER SETUP PROGRAM 

This program reads the AEGL line listing tape and creates a file of 
selected line parameters. The two-step JCL routine runs the FORTRAN 
program to read and select the line parameters and then calls for an all 
core sort using the IBM system sort routine to put the lines in increasing 
order of the lowest frequency value at which they have significant absorp- 
tion. A sample of the line parameter output after sorting is also in- 
cluded. As written, the program selects N 2 O lines from 2550 - 2560 cm -1 
for a ground-based observer looking straight up at the sun. The line 
parameter listing example is for a different case. 
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